Dynamics and Synchrony from Oscillatory Data via Dimension Reduction 
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Complex, oscillatory data arises from a large variety of biological, physical, and social systems. 
However, the inherent oscillation and ubiquitous noise pose great challenges to current methodology 
such as linear and nonlinear time series analysis. We exploit the state of the art technology in 
pattern recognition and specifically, dimensionality reduction techniques, and propose to rebuild 
the dynamics accurately on the "cycle" scale. This is achieved by deriving a compact representation 
of the cycles through global optimization, which effectively preserves the topology of the cycles that 
are embedded in a high dimensional Euclidian space. Our approach demonstrates a clear success in 
capturing the intrinsic dynamics and the subtle synchrony pattern from uni/bivariate oscillatory data 
over traditional methods. Application to the human locomotion data reveals important dynamical 
information which allows for a clinically promising discrimination between healthy subjects and 
those with neural pathology. Our results also provide fundamental implications for understanding 
the neuromuscular control of human walking. 

PACS numbers: 05.45.Tp, 05.45.Xt 
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Time series exhibiting complex oscillatory dynamics 
are widely observed in the real world: from finance to 
biomedical systems [l| . In this Letter we focus on an im- 
portant class of data with a stable oscillation frequency 
but irregular waveform fluctuations, also called pseu- 
doperiodic time series. Such data arise from broad ap- 
plication domains and have gained particular interest in 
recent years, with examples ranging from human ECG 



and gait data 
data in patients 




white blood-cell count and tremor 
^ ^ , epidemic dynamics ^6^ , light inten- 
sity of laser , sun spot numbers and global temperature 
variation [s, 9]. Accurate extraction and characterization 
of the dynamics of these time series is a general problem, 
which holds the key to understanding the inner work- 
ings of many important physical and biological systems. 
To this end, traditional methods rely primarily on linear 
spectral analysis or computing nonlinear characteristics, 
and recent attempts include producing pseudoperiodic 
surrogate data [lC| , or to establish a nonlinear transform 
from cycles in time domain to nodes in a dual complex 
network domain [ll| . 

However, there are still no generic, systematic and ro- 
bust approaches to handle such oscillatory time series. 
The Fourier transform, though widely applied to oscil- 
latory data, is inherently a linear method and cannot 
account for the nonlinear nature of the data, if any. For 
nonlinear approaches, the cyclic trend overlying the sig- 
nal, together with noise from unknown sources, can of- 
ten mask the intrinsic dynamics 12, and pose great 
challenges to the nonlinear time series analysis/modeling 
techniques. For example, the most popular Poincare sec- 
tion method, which reduces the flow data to intersec- 
tions of the trajectories with a lower-dimensional sub- 
space, may produce highly corrupted results under noisy 
environment, in that the intersection points can no longer 



preserve the original dynamics in the presence of noise. 
Things get even worse when the cycles possess complex 
noncoherent waveforms. In this paper we propose a 
novel, generic approach that can effectively capture the 
dynamics of oscillatory time series. By mapping the cy- 
cles in the time series to a multi-dimensional Euclidian 
space, we seek a low-dimensional representation which 
topologically preserves the important proximity relation 
among cycles, through efficient global optimization. It 
is the first time that advances in dimension reduction 
are explored in reconstructing the dynamics of complex 
oscillatory data. Our approach utilizes the richer infor- 
mation of pairwise cycle correlation, therefore it not only 
excavates the inherent dynamics obscured by the cyclic 
trend, but also offers an extra robustness to noise due 
to the global nature of the method. Exploration of our 
approach to probing the synchrony between bivariate os- 
cillatory time series is shown to be able to reveal the sub- 
tle synchrony for which traditional approaches will fail. 
We applied the proposed method to the human locomo- 
tion data, and are capable of discriminating between the 
healthy subjects and those with neuropathology reliably. 
Based on the results we are able to make the important 
biological inference that the human walking is not criti- 
cally dependent on the peripheral neural feedback. 



THEORY 

Reconstructing Dynamics Underlying Cyclic Trend 
by Spectral Clustering 

We demonstrate how intrinsic dynamics of psuedope- 
riodic data can be extracted using Laplacian Eigen- 
maps |l4j , a nonlinear dimension reduction method based 
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FIG. 1: (a) Time series form X-component of the chaotic 
Rossler system, which is divided into cycles by local minima 
(circles), (b) Laplacian embedding c(i) extracted on the cycle 
scale for the time series shown in (a), and (c) the return plots 
of c(i) derived by two dimension reduction methods, i.e., the 
Laplacian Eigenmaps and Multidimentional scaling. 



FIG. 2: Return plots for the time series from (a) Poincare 
section points P{i), which is obtained by collecting the local 
minima Xmin , and (b) the extracted c(i) . The Rossler system 
here is corrupted by 5% dynamical noise and 30% measure- 
ment noise. It is obvious that the return plots from c(i) takes 
on a clearer form of quadratic function, indicating that the 
original nonlinear dynamics is sufficiently kept. 



on spectral graph theory. We illustrate with the X- 
component of the chaotic Rossler system. Motivated by 
the fact that such data usually exhibit a highly redundant 
patterns in the form of repeated cycles, we first partition 
the time series into individual cycles C^'s {i — 1, k) by 
local minimums as shown in Fig. [T] (a) . Each cycle is then 
mapped to a high dimensional vector Xi whose dimension 
equals the number of points in that cycle. Our goal is 
to compute a set of new, low-dimensional (in the sim- 
plest case, ID) representation, or embedding y^'s, such 
that the proximity relations among x^'s are maximally 
preserved in their low-dimensional counterparts j/i's and 
that the redundancy in the cycles are removed. 

To achieve this, a weighted graph ^ is constructed, 
where each node corresponds to a cycle Xi and edges are 
assigned between all pairs of nodes. The graph can ei- 
ther be fully connected or only bind those vertices within 
the fc-nearest-neighbors and we aims at extracting the 
dynamics from the topology of the graph. We use Wij 
to denote the similarity between cycle Xi and Xj, which 
can be chosen conveniently as the correlation coefficient 
Pij = Cov{xi, Xj)/{axiCrxj) (in case Xi and Xj differ in 
length, shift the shorter vector on the longer one un- 
til Pij maximizes). Then, the low-dimensional represen- 
tation 2/i's can be cast as the solution of the following 
optimization problem, min^^ Wij — J/j|P, which 
penalizes those mappings where nearby points 
relocated far apart in the space of j/i's. In case of uni- 
variate 2/i's, the objective can be written as y^Ly, where 
y = [2/1,2/2, •■•,2/fc]^, L = D - W is the graph Lapla- 
cian, and D is the diagonal degree matrix such that 
Da — J2j ^ji- Hence the name Laplacian Eiganmap. 
To prevent arbitrary scaling in the solution, one can en- 
force the constraint y^Dy = 1. 



The above constrained minimization is solved by the 
generalized eigenvalue problem Ly^ = XiDy^, where A^'s 
(« = 1, 2, fc) are eigenvalues sorted in an ascending or- 
der, and y/s are the corresponding eigenvectors. It can 
be shown that the minimum eigenvalue Ai is zero, cor- 
responding to an eigenvector (y^) of all I's. Therefore 
it is degenerate and the optimal solution is actually pro- 
vided by Y t, the eigenvector of the second smallest eigen- 
value 1J|. As is shown in Fig. [1] the eigenvector y2(i) 
provides a unique, reduced representation of the origi- 
nal time series by capturing the dynamics of the oscilla- 
tory data on the cycle scale. We use a general notation 
c(i) (c for cycle) for such reduced series representation 
(i.e., other dimension reduction schemes can also be ap- 
plied, and we denote the reduced series all as c(i)). More 
generally, to obtain an m-dimensional solution for y^'s, 
one can simply extract the m trailing eigenvectors y/s 
(i = 2,3,...,m-Hl). 



The Significance of Extracting Dynamics on Cycle 
Scale 



The comparison between Laplacian embedding c{i) 
and the Poincare section points P{i) (obtained by collect- 
ing the local minima) indicates that they are dynamically 
identical by sharing the same chaotic invariants, i.e., Lya- 
punov exponent and correlation dimension. However, in 
case of significant noise (see Fig. [2]), our approach can 
capture more of the deterministic structure than does 
P{i) and is therefore more robust. This is not surpris- 
ing, since acquisition of c(i) is based on an optimization 
process that preserves the proximity relation between all 
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pairs of cycles simultaneously, while P{i) is obtained by 
treating each cycle independent of each other. It is rea- 
sonable to expect the former to excavate more useful, 
global patterns than the latter. 

Many other forms of dimension reduction can also be 
applied here, which produce similar results. For example, 
the Multidimensional Scaling method that preserves the 
pairwise distance yields almost the same outputs as the 
Laplacian Eigenmaps, see Fig. [Ijc). Actually, the appli- 
cability of dimension reduction techniques in extracting 
the dynamics shall be generally justifiable, considering 
the intrinsically low correlation dimensions of most real 
world pseudoperiodic data. The corresponding trajecto- 
ries in phase space tend to have similar orientations for 
nearby cycles (see, e.g.. Fig. SJ. Such redundancy can be 
effectively removed through dimension reduction, leav- 
ing only the useful degrees of freedom for manifestation 
of the dynamics of interest. 



5 



(b) 














r 



P,(i) c,(i) 

FIG. 3: Correlation between X and Y component revealed by 
(a) Poincare section points P{i) and (b) the extracted c(i). 
The Rossler system is corrupted by 5% dynamical noise and 
30% measurement noise. 



Detecting Synchrony from Oscillatory Data 

A topic which is of great interest to oscillatory data 
analysis in recent years is to detect the degree of syn- 
chronization between self-sustained oscillators. For ex- 
ample, the peakness of the phase difference distribution 
and the consistency of mutual nearest neighbors are used 
to characterize the phase synchrony Il7l | and the dy- 
namical interdependence [l8|, respectively. Here we are 
especially interested in the case where the phase relation 
between two oscillator is evident, but is hard to estimate 
due to non-phase-coherence and noise, or is not sensi- 
tive to degree of changes of synchrony. For example, the 
blood pressure interacted with heartbeat, where each cy- 
cle of blood pressure variation correspond to one heart- 
beat. The phase index in this case may not be quite infor- 
mative of the synchrony strength. On the other hand, the 
presence of noise in most real world data will inevitably 
degrade the dependency measures. 

Unlike traditional methods, we propose to quantify the 
synchrony the noisy, noncoherent pseudoperiodic data 
on theAr cycle scale through their Laplacian embedding 
c{i). We first test with the X and Y components of 
the noisy Rossler system. Although these two compo- 
nents are non-separable and cannot be strictly defined 
as being synchronized, they are actually "in phase" and 
therefore suitable in our context. The two time series 
are first segmented into cycles by the local maximums 
of one of the data (i.e., the zth cycle in both time series 
share the same time span). We find that the Poincare 
section points picked up from the two time series can 
hardly capture the interdependency due to noise, in that 
the scatter plot between the Poincare section points of 
X and Y leads to a group of random points (Fig. [3] (a)). 
Comparatively, the extracted Laplacian embedding c(i)'s 
extracted from X and Y by our approach successfully re- 
veal the synchrony pattern by demonstrating an uprising 
trend in the corresponding scatter plot (Fig. [3] (b)). 



Reconstructing the dynamics on the cycle scale not 
only enhances the robustness to noise, but bring new vi- 
tality to a number of nonlinear methods which are other- 
wise not suitable for oscillatory data due to the inherent 
periodicity, such as detrended fluctuation analysis, sur- 
rogate data method, recurrence plot, causality, entropy 
and so on. For example, the detrended fluctuation anal- 
ysis, which computes the variance of the detrended data 
at different scales, will produce similar variance curves 
(i.e., uprising at small scale and saturate at the periodic- 
ity of the signal) for all pseudoperiodic data despite their 
distinct dynamics. This is because the periodic trend 
dominating in the data will conceal the underlying, sub- 
tle dynamics. Thus extracting the dynamics on the cycle 
scale proves to be crucial for subsequent analysis. 



APPLICATION TO HUMAN LOCOMOTION 
DATA 

Characterizing the Locomotion Dynamics 

Now we apply our approaches to human gait data col- 
lected from electrogoniometers at the knee and ankle 
joints that measure their sagittal plane kinematics dur- 
ing overground walking. Human locomotion is a highly 
complex, rhythmic process that involves multilevel con- 
trol of central nerve system and feedback from various 
peripheral sensors. Here we consider two groups of sub- 
jects from healthy controls (CO) and diabetic patients 
(NP, with significant peripheral neuropathy), each with 
10 subjects 3]. 

Typically, the human gait time series (Fig. [5] (a) and 
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FIG. 4: Time series (left panel) and the corresponding phase 
space reconstructions (right panel) for knee ((a)-(b)) and an- 
kle kinematics data ((c)-(d)). The two time series are typ- 
ically non-phase-coherent by demonstrating multi-oscillation 
within each cycle. This is also evident from the multi-center 
rotations of the attractor in phase space (see (b) and (d)). 
The two signals are divided by their respective local maxima 
into consecutive cycles. 



large time span are still statistically correlated). In com- 
parison, the spectrum of the diabetic patients are mostly 
flat resembling white noise processes (/3 = 0.37 ± 0.16), 
which means that the strides at different times are mostly 
uncorrelated. This difference, however, has not been 
found with either the stride interval (SI) series (Fig. [S] 
(c)-(d)) or the raw data (Fig. [5](e)-(f)). This is because, 
the periodicity in the raw data can obscure the under- 
lying fluctuations and that a linear Fourier transform is 
not capable of characterizing nor distinguishing the non- 
linear dynamics in the data; on the other hand, although 
SI removes the cyclic trend and has been widely used to 
quantify pathological locomotion ^iB'l , it loses much of the 
dynamical information by only recording the duration of 
cycles. In comparison, the transform c(i) successfully re- 
moves the periodic trend while preserving the dynamical 
fluctuation within each cycle, in particular, the specific 
patterns of the four phases within a pace. Therefore it 
keeps more valuable information about the neuromuscu- 
lar control of human walking. 
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FIG. 5: Power spectrum density (PSD) for typical healthy 
subject (the left panel) and the diabetic patient (the right 
panel). The top, middle and bottom rows are PSDs for the 
extracted c(i), stride interval series, and the original ankle 
data, respectively. The frequency here is normalized by the 
sampling rate. The PSDs for the stride interval series and 
the original data show no significant difference between the 
healthy subject and the diabetic patient. 



(c)) exhibit a relatively fixed frequency with the stride 
interval demonstrating little variability, which is suited 
for dimension reduction. We first extract the dynam- 
ics of the ankle on its cycle scale by Laplacian em- 
bedding. Figure [5] (a)-(b) display the power spectrum 
density (PSD) of the extracted c{i) typically observed 
for the two groups. We find that most CO subjects 
demonstrate broad band spectrums that scale as l/f^, 
P = 0.76{mean) ± 0.23(std), indicating the presence of 
long range correlation (i.e., the strides separated by a 



Synchrony Detection between Knee and Ankle 
Movement 



The human walking involves the coordination of two 
major joints, i.e., the knee and the ankle. So we are also 
interested in examining the synchrony between the knee 
and ankle data, which may provide further insights in 
understanding the neuralcontrol of walking. The knee 
and ankle movements during continuous walking are ob- 
viously in phase due to the physical connection of the two 
joints. The strength of coupling, however, can hardly 
be recognized from the phase due to the noncoherence 
[lot (see Fig. S]). Also, noise tends to destroy the local 
structure in phase space and thus hampers the dynami- 
cal dependency measures [2^. To circumvent these diffi- 
culties, we compare the dynamics of the two time series 
by using their Laplacian Embedding c(i)'s. Note that 
each time series can be segmented by either its own lo- 
cal maximums, or those of its partner series (markers in 
Fig. [5]) . Therefore we will segment each time series twice 
and compute the averaged correlation coefficients pij^s 

between Cankle and Cknee- 

Then we examine the interrelation between the Lapla- 
cian embeddings from the ankle and the knee data. For 
healthy subjects, the scatter plots demonstrate a signif- 
icant uprising trend (Fig. [U (a), p — 0.68 ± 0.19), in- 
dicating that the knee and ankle movements are highly 
synchronized; while for diabetics patients, the synchrony 
almost vanishes (Fig. [6] (b), p = 0.26 ± 0.18). Again, 
the discrimination between CO and NP groups is un- 
available by stride interval series, which always exhibits 
a strong correlation between the joints (Fig. [S] (c)-(d)), 
corresponding to high degree of phase synchronization. 
Note that our approach avoids the difficulty of defining 
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FIG. 6: Synchrony between knee and ankle locomotion re- 
vealed by Laplacian embedding c(i) ((a) and (b)) and stride 
interval series ((c) and (d)), with the left and right columns 
for an healthy subject and a diabetic patient, respectively. 
As can be seen in (b), the stride interval series cannot distin- 
guish the healthy from the diabetics, as the scatter plots of 
the stride interval series show similar uprising trends for both 
of them. 



the phase. On the other hand, it is superior to directly 
computing the correlation between pairs of equal-time cy- 
cles picked respectively from the two time series, which is 
sensitive to noise and nonstationarity. Finally, we point 
out that a more comprehensive description of synchrony 
can be achieved by examining more Laplacian eigenvec- 
tors. In the current case the single eigenvector already 
encodes the primary variability and is thus sufficient for 
the discriminative task. 

Note that we have observed a lack of long range correla- 
tion in the ankle kinematics of the NP group (Fig. [5](b)). 
This suggests the alteration of the locomotor pattern in 
the lower limbs, due to the nerve deterioration in feet 
and toes typical of the diabetics. Despite the impaired 
peripheral feedback caused by the dying nerves, the knee 
kinematics of the NP group are still found to demonstrate 
a stable long range correlation indistinguishable from the 
CO group. This "mismatch" between the ankle/knee dy- 
namics for diabetics can be reasonably explained by the 
weak synchrony between these two joints, as is shown in 
Fig. [S] (b). From these findings, we can see that the im- 
paired neuralfeedback from the feet influences only the 
lower limb locomotion while not that of the knees. We 
therefore conclude that the human walking is not criti- 
cally dependent on the feedback from peripheral nervous 
system. It should be noted that although Gates et.al. 
2l\ checked the same data, they do not consider the in- 
teraction among the knee and angle locomotion which 
may be important in characterizing understanding neu- 
ralcontrol of human walking. Furthermore the authors 
relied solely on the extracted stride interval series, which 
have been shown to contain limited dynamical informa- 



tion and therefore may be insufficient for the analysis and 
evaluation of the human walking. 



CONCLUSION 

To conclude, we have for the first time circumvented 
the problem of rebuilding the complex dynamics in oscil- 
latory data from the viewpoint of dimension reduction. 
We use a global optimization procedure to compute the 
inherent, low-dimensional representation that maximally 
preserves the topology of the cycles when embedded in 
the multi-dimensional Euclidian space. Our approach 
has been shown to be very robust to noise, and is capa- 
ble of extracting the underlying dynamics and identifying 
the subtle synchrony pattern which usually degrade tra- 
ditional linear or nonlinear methods. Application to hu- 
man gait data provides clinically promising information 
in discriminating the healthy from the neuropathologi- 
cal patients, and further enables us to make fundamental 
inference on the neuromuscular control mechanism of hu- 
man walking. 

Our approach may be of great relevance, and is ex- 
pected to provide more accurate and powerful diagnos- 
tics, to the general biological and biomedical fields where 
complex oscillatory data abound. For example, we have 
applied the proposed method to the phonation data from 
the Parkinson patients, and find a significant lack of long 
range correlation among the cycles in the signal com- 
pared to the healthy people. Our approach can therefore 
serve as an alternative tool in early Parkinson disease 
(PD) detection, where there are no blood or laboratory 
tests currently that can help in diagnosing PD. Other 
potential applications which will be involved in the fu- 
ture works include the discrimination between Parkinson 
and essential tremor, and the evaluation of the interac- 
tion between the blood pressure variation with the heart 
beat. 
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